/*
* CIEDE2000.cpp
* Part of http://github.com/gfiumara/CIEDE2000 by Gregory Fiumara.
* See LICENSE for details.
*/
#include "stdafx.h"
#include <cmath>

#include "CIEDE2000.h"

/*******************************************************************************
* Conversions.
******************************************************************************/

double CIEDE2000::deg2Rad(const double deg)
{
	return (deg * (M_PI / 180.0));
}

double CIEDE2000::rad2Deg(const double rad)
{
	return ((180.0 / M_PI) * rad);
}

void CIEDE2000::RGB2LAB(const RGBQUAD *pRGB, LAB *pLAB)
{
	int LabTable[1024];
	for (int i = 0; i < 1024; i++)
	{
		if (i > 9)
			LabTable[i] = (int)(pow((float)i / 1020, 1.0f / 3) * (1 << 10) + 0.5);
		else
			LabTable[i] = (int)((29 * 29.0 * i / (6 * 6 * 3 * 1020) + 4.0 / 29) * (1 << 10) + 0.5);
	}
	const static int big_shift = 18;
	const static int HalfShiftValue = 512;
	const static int shift = 10;
	const static int offset = 128 << shift;
	const static int ScaleLC = (16 * (1 << shift));
	const static int ScaleLT = 116;
	const static int ScaleY = 903;
	const static int para1 = 500;
	const static int para2 = 200;
	const static int ThPara = 9;
	long long X = (pRGB->rgbBlue * 199049 + pRGB->rgbGreen * 394494 + pRGB->rgbRed * 455033 + 524288) >> (big_shift);
	long long Y = (pRGB->rgbBlue * 75675 + pRGB->rgbGreen * 749900 + pRGB->rgbRed * 223002 + 524288) >> (big_shift);
	long long Z = (pRGB->rgbBlue * 915161 + pRGB->rgbGreen * 114795 + pRGB->rgbRed * 18621 + 524288) >> (big_shift);

	pLAB->l = Y > ThPara ? ((ScaleLT * LabTable[Y] - ScaleLC + HalfShiftValue) >> shift) : ((ScaleY* Y) >> shift);
	pLAB->a = ((para1*(LabTable[X] - LabTable[Y]) + HalfShiftValue + offset) >> shift) - 128;
	pLAB->b = ((para2*(LabTable[Y] - LabTable[Z]) + HalfShiftValue + offset) >> shift) - 128;
	/*
	double x = 0.412453*pRGB->rgbRed + 0.357580*pRGB->rgbGreen + 0.180423*pRGB->rgbBlue,
		y = 0.212671*pRGB->rgbRed + 0.715160*pRGB->rgbGreen + 0.072169*pRGB->rgbBlue,
		z = 0.019334*pRGB->rgbRed + 0.119193*pRGB->rgbGreen + 0.950227*pRGB->rgbBlue;
	x = x / (255.0*0.950456);
	y = y / 255.0;
	z = z / (255.0*1.088754);
	double fx, fy, fz, l, a, b;
	if (y>0.008856)
	{
		fy = pow(y, 1.0 / 3.0);
		l = 116.0*fy - 16.0;
	}
	else
	{
		fy = 7.787*y + 16.0 / 116.0;
		l = 903.3*y;
	}
	if (x>0.008856)
	{
		fx = pow(x, 1.0 / 3.0);
	}
	else
	{
		fx = 7.787*x + 16.0 / 116.0;
	}
	if (z>0.008856)
	{
		fz = pow(z, 1.0 / 3.0);
	}
	else
	{
		fz = 7.787*z + 16.0 / 116.0;
	}
	a = 500.0*(fx - fy);
	b = 200.0*(fy - fz);
	pLAB->l = l;	//0-100
	pLAB->a = a;	//-128-127
	pLAB->b = b;	//-128-127
	*/
}

double CIEDE2000::CIEDE2000(const LAB &lab1, const LAB &lab2)
{
	/*
	* "For these and all other numerical/graphical ??delta E00 values
	* reported in this article, we set the parametric weighting factors
	* to unity(i.e., k_L = k_C = k_H = 1.0)." (Page 27).
	*/
	const double k_L = 1.0, k_C = 1.0, k_H = 1.0;
	const double deg360InRad = CIEDE2000::deg2Rad(360.0);
	const double deg180InRad = CIEDE2000::deg2Rad(180.0);
	const double pow25To7 = 6103515625.0; /* pow(25, 7) */

	/*
	* Step 1
	*/
	/* Equation 2 */
	double C1 = sqrt((lab1.a * lab1.a) + (lab1.b * lab1.b));
	double C2 = sqrt((lab2.a * lab2.a) + (lab2.b * lab2.b));
	/* Equation 3 */
	double barC = (C1 + C2) / 2.0;
	/* Equation 4 */
	double G = 0.5 * (1 - sqrt(pow(barC, 7) / (pow(barC, 7) + pow25To7)));
	/* Equation 5 */
	double a1Prime = (1.0 + G) * lab1.a;
	double a2Prime = (1.0 + G) * lab2.a;
	/* Equation 6 */
	double CPrime1 = sqrt((a1Prime * a1Prime) + (lab1.b * lab1.b));
	double CPrime2 = sqrt((a2Prime * a2Prime) + (lab2.b * lab2.b));
	/* Equation 7 */
	double hPrime1;
	if (lab1.b == 0 && a1Prime == 0)
		hPrime1 = 0.0;
	else {
		hPrime1 = atan2(lab1.b, a1Prime);
		/*
		* This must be converted to a hue angle in degrees between 0
		* and 360 by addition of 2?? to negative hue angles.
		*/
		if (hPrime1 < 0)
			hPrime1 += deg360InRad;
	}
	double hPrime2;
	if (lab2.b == 0 && a2Prime == 0)
		hPrime2 = 0.0;
	else {
		hPrime2 = atan2(lab2.b, a2Prime);
		/*
		* This must be converted to a hue angle in degrees between 0
		* and 360 by addition of 2?? to negative hue angles.
		*/
		if (hPrime2 < 0)
			hPrime2 += deg360InRad;
	}

	/*
	* Step 2
	*/
	/* Equation 8 */
	double deltaLPrime = lab2.l - lab1.l;
	/* Equation 9 */
	double deltaCPrime = CPrime2 - CPrime1;
	/* Equation 10 */
	double deltahPrime;
	double CPrimeProduct = CPrime1 * CPrime2;
	if (CPrimeProduct == 0)
		deltahPrime = 0;
	else {
		/* Avoid the fabs() call */
		deltahPrime = hPrime2 - hPrime1;
		if (deltahPrime < -deg180InRad)
			deltahPrime += deg360InRad;
		else if (deltahPrime > deg180InRad)
			deltahPrime -= deg360InRad;
	}
	/* Equation 11 */
	double deltaHPrime = 2.0 * sqrt(CPrimeProduct) *
		sin(deltahPrime / 2.0);

	/*
	* Step 3
	*/
	/* Equation 12 */
	double barLPrime = (lab1.l + lab2.l) / 2.0;
	/* Equation 13 */
	double barCPrime = (CPrime1 + CPrime2) / 2.0;
	/* Equation 14 */
	double barhPrime, hPrimeSum = hPrime1 + hPrime2;
	if (CPrime1 * CPrime2 == 0) {
		barhPrime = hPrimeSum;
	}
	else {
		if (fabs(hPrime1 - hPrime2) <= deg180InRad)
			barhPrime = hPrimeSum / 2.0;
		else {
			if (hPrimeSum < deg360InRad)
				barhPrime = (hPrimeSum + deg360InRad) / 2.0;
			else
				barhPrime = (hPrimeSum - deg360InRad) / 2.0;
		}
	}
	/* Equation 15 */
	double T = 1.0 - (0.17 * cos(barhPrime - CIEDE2000::deg2Rad(30.0))) +
		(0.24 * cos(2.0 * barhPrime)) +
		(0.32 * cos((3.0 * barhPrime) + CIEDE2000::deg2Rad(6.0))) -
		(0.20 * cos((4.0 * barhPrime) - CIEDE2000::deg2Rad(63.0)));
	/* Equation 16 */
	double deltaTheta = CIEDE2000::deg2Rad(30.0) *
		exp(-pow((barhPrime - deg2Rad(275.0)) / deg2Rad(25.0), 2.0));
	/* Equation 17 */
	double R_C = 2.0 * sqrt(pow(barCPrime, 7.0) /
		(pow(barCPrime, 7.0) + pow25To7));
	/* Equation 18 */
	double S_L = 1 + ((0.015 * pow(barLPrime - 50.0, 2.0)) /
		sqrt(20 + pow(barLPrime - 50.0, 2.0)));
	/* Equation 19 */
	double S_C = 1 + (0.045 * barCPrime);
	/* Equation 20 */
	double S_H = 1 + (0.015 * barCPrime * T);
	/* Equation 21 */
	double R_T = (-sin(2.0 * deltaTheta)) * R_C;

	/* Equation 22 */
	double deltaE = sqrt(
		pow(deltaLPrime / (k_L * S_L), 2.0) +
		pow(deltaCPrime / (k_C * S_C), 2.0) +
		pow(deltaHPrime / (k_H * S_H), 2.0) +
		(R_T * (deltaCPrime / (k_C * S_C)) * (deltaHPrime / (k_H * S_H))));

	return (deltaE);
}
